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Abstract 

The relative computational effort among the spatially five-point numerical flux functions of 
Harten, van Leer, and Osher and Chakravarthy is explored. These three methods typify the 
design principles most often used in constructing higher than first-order upwind total variation 
diminishing (TVD) schemes. For the scalar case the difference in operation count between any 
two algorithms may be very small and yet the operation count for their system counterparts 
might be vastly different. The situation occurs even though one starts with two different yet 
equivalent representations for the scalar case. 


I. Introduction 

Various high-resolution upwind difference schemes have been developed in recent years [1- 
9j. Of special interest are the extension of these methods to systems of nonlinear hyperbolic 
conservation laws in one and higher dimensions. Typical methods for systems are MUSCL 
[1,2], Osher [6], Harten [3], Roe [5], and van Leer [1], Here van Leer method means using 
his methodology [1] to extend the first-order Roe’s scheme [5] to higher than first-order and 
MUSCL means van Leer’s pioneer work on extending the first-order Godunov scheme to higher- 
order. The MUSCL and Osher type of approach require considerably more computation and 
programming effort than with the latter three, especially for problems in multidimensions and in 
curvilinear coordinates. Harten and Roe’s second-order schemes share many common features. 
Considering only the latter three, the way second-order accuracy is achieved and the form of 
limiters play an important role. For the scalar case, the difference in operation count between 
any two algorithms may be very small and yet the operation count for their system counterparts 
might be vastly different. This situation occurrs even though one starts with two different yet 
equivalent representations for the scalar case. An example will be given in section V. 

Recently, Osher and Chakravarthy [10] extended their earlier work on five-point fully upwind 
second-order total variation diminishing (TVD) schemes to five-point second-order or third-order 
upwind biased TVD schemes. They also derived formulas for more than five points and higher 
than third-order cases. Only the second and third-order are discussed here. They claimed that 
their five-point second-order or third-order upwind biased scheme requires hardly anymore work 
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than their five-point fully upwind second-order scheme. It was the investigation of their work 
which prompted the work of this note. As one can see later for system cases, their five-point 
upwind biased schemes require more operations than their fully upwind second-order scheme. 
Furthermore, their five-point fully upwind or upwind biased schemes require substantially more 
operations than the Harten or van Leer method. For the purpose of this paper, order of accuracy 
of a TVD scheme refers to the numerical solution that is away from extreme points. 

Specifically, it is the objective of this note to discuss the relative computational effort among 
the spatially five-point TVD schemes of Harten, van Leer and Osher and Chakravarthy [3,7,10] 
for system cases via the Huang [11] and Roe’s [12] generalization for systems. It is emphasized 
here that, the scalar schemes of Harten, van Leer and Osher and Chakravarthy are used. Then 
extensions to systems are employed by the same method (i.e. Huang or Roe), not by their 
original Riemann solver [1,13] for system cases. Huang and Roe’s generalization to system cases 
is based on a local linearization which defines at each point a “local” system of characteristic 
fields, and then applies the scalar scheme to each of the scalar characteristic equations. Because 
of the freezing of coefficients, one can transform the local characteristic field back to the original 
field variables, and arrive at an algorithm for the nonlinear hyperbolic system case. Numerical 
experiments with these schemes on one and two-dimensional steady-state calculations [14-20] 
show that they produce similar shock resolution. In certain instances, the third-order scheme 
gives a slightly better result than the second-order method [20], In general, the improvement 
from second-order to third-order is far less pronounced than from first-order to second-order 
TVD schemes [19], This is due partly to the fact that all TVD schemes of higher than first- 
order reduce to first-order at points of extrema. Knowing their performance, it is relevant to 
compare their computational effort. 

The comparison is only on the numerical flux function. No specific time discretization is 
involved. For simplicity in presentation, the forward Euler time differencing is assumed. One 
should notice that under the forward Euler time differencing only Harten’s method is second- 
order in time and space; the other two are first-order in time and second-order in space. As 
a side remark, TVD schemes are not necessarily upwind [21-23], but upwind TVD schemes, in 
general, give sharper shocks than symmetric TVD schemes. The advantage of symmetric TVD 
schemes is that they are easier to implement and provide a more natural way of extending the 
scheme to two and three-dimensional problems. One can easily modify Harten’s method to be 
non-upwind. In this note, only the aforementioned three upwind methods are considered. 


II. Harten TVD Scheme 

Consider a scalar hyperbolic conservation law 

( 2 . 1 ) 

where / is the flux function and a(u) = df/du is the characteristic speed. Let u” be the 
numerical solution of (2.1) at x — jAx and t = nAt, with Ax the spatial mesh size and At the 
time step. Consider a general five-point explicit difference scheme in conservation form 


^ = 0 
dt dx 
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( 2 . 2 ) 


.n+1 

l 3 


A- *?_*)• 


where /i" ±jL = u",u" ±1 , u ^± 2 )> anc ^ A = At/Ax. Here /i, commonly called a numerical 

flux function, is required to be consistent with the conservation law in the following sense 

h(u,u,u,u) = f(u). (2.3) 

The numerical flux of the second-order explicit TVD scheme of Harten [3] can be expressed 


7h 

h , + i . 


\\fj + fj+i - ip(a 3+ ±)A i+ xu], 


where f, = f{u 3 ) + gj, a j+ i = a j+ i + 7 y+ i, and A j+ iu = u y+1 - Uj , with 

xl>(z) = \z\ 




( [fj +1 /j)/A J+ l U Ay + iU ^ 0 


1 «K) 




u — 0 


(2.4a) 

(2.4b) 

(2.4c) 


and 


9j = s ■ max [°. m in( cr j + i I A i+ 1 u|, S • a 3 _ i A J _ i u)] 
S — sgn(A J+ i u) 


'liH 


(#+i - 3 y)/A y+ iu 
0 




iu ^ 0 
u = 0. 


Here Cy + i = &{ a j+ a) and can be expressed as 


<?( z ) = “ Az2 l ^ 0 


(2.4d) 


(2.4e) 


(2-4*) 


for time accurate calculations. With this choice of (t(z), the scheme (2.2) together with (2.4) 
is second-order accurate in both time and space (except at points of extrema); see [3,4]. Also, 
with this choice of the g function, the scheme will automatically switch itself to first-order at 
points of extrema. This is one of the vehicles to prevent spurious oscillations near a shock. 

Define a “minmod” function of a list of arguments to be the smallest number in absolute 
value if the list of arguments are of same sign, but to be zero if any argument is of the opposite 
sign. One can rewrite (2.4d) as gj = minmod(ffy + i Aj + iu,aj_x Aj_xu). It is a form of the 
so called “limiter” for the control of unwanted oscillations in numerical schemes. Other more 
general forms for g can be obtained following the line of argument of Roe [5]. As a side remark, 
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if one sets all the fly’s to zero in (2.4), the resulting scheme becomes Roe’s first-order upwind 
method. 

With the choice of the ip function in(2.4b), the corresponding first-order upwind scheme is not 
consistent with an entropy inequality, and the scheme might converge to a nonphysical solution. 
A slight modification of the numerical viscosity term 

^ = { (J 2 + e 2 )/2e jzj < e 

can remedy the entropy violating problem [3,4], where e > 0 is a parameter. A formula for 
e can be found in reference [24]. The ip{z) in (2.5) is a continuously differentiable positive 
approximation to |z| in (2.4b). 

For simplicity of presentation, assume ip(z) has the form (2.4b). The use of an entropy 
satisfying ip(z) is only a minor modification of what is about to be presented. Furthermore, 
since the main concern of this note is the spatial discretization, only the forward Euler time 
differencing is considered. Other explicit, implicit, or predictor and corrector types of time 
differencing do not alter the results. 

Using the identity |ay + i| = sgn(Sy + i)a J+ i and ay + i = (/y+i - /y+i)/(«j+i - uy), Marten’s 
scheme can be rewritten as 

u " + ‘ = «? - J [' - ‘S-'Pi+j)] WVi - !j ) - | [l + s8»(S“.j)] (/? ' I?- 1). (2.6) 

Here, the numerical flux function becomes 

kf+i = \ [* ~ sgn(a i+ i)] ifJ+i - U) + U- (2-7) 

Note that with the choice of g 3 given in (2.4d), sgn(a" +i ) = sgn(a" +i ). Thus, one does not 
have to calculate 7 J+ i. The h^ +L can also be expressed as 

h?+i = \ [i - sgn(a i+ i)] (f J+ 1 - fj) + fj. (2.8) 

With the use of the numerical flux (2.7) as an alternative, the g function of (2.4d) can be 
defined in a slightly different form. 

9i = minmod (^ + iA J + i/,^.iA 3 .i/) (2.9a) 


with <Ty + i = a(a^i) and 

°( z ) = \ [»gn(*) - A*]. (2.9b) 

Here the identities sgn(A J + i/)sgn(a 3 + i) = sgn(u J + 1 - uy) and |A J+ i/| = |oy + i ||uy +1 - uy| 
are used. The above limiter g : of (2.9) can be considered as a flux limiter since the flux / is 
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limited. For scalar problems, equations (2.4a), (2.7) and (2.8) are equivalent. The difference in 
operation count are very minor. Equations (2.4d,f) are preferred over (2.9) for its straightforward 
extension to system cases because u appears rather than /. Furthermore, in the extension to 
nonlinear systems via Huang’s and Roe’s generalizations, (2.4a) results in a smaller operation 
count than (2.7) or (2.8). 


III. van Leer Differencing 

In references [7,19,25] van Leer and his collaborators adopted the original MUSCL [1] type 
of differencing to obtain a second and third-order flux vector splitting method. Very impressive 
results were obtained for steady-state airfoil calculations. Their result indicated that one does 
not have to resort to a Godunov type of Riemann solver in order to obtain good shock resolution. 
In this paper the application of the MUSCL type differencing to Roe’s first-order scheme is 
proposed. 

Consider a three point explicit difference scheme in conservation form 


= »" - - *(«?-..«”)] ( 3 -l) 

where = A(u”, u" +1 ), is a first-order numerical flux. Here 1 without a “ “ ” is used 

to distinguish a first-order numerical flux from the higher-order numerical flux + 

Roe’s first-order scheme (without entropy fixed) can be written as 
^j + 3 = 2 I/*”*' fc +1 ~ l a i+3l^i+5 U ]’ 


with defined in (2.4c). 

The MUSCL approach to obtain a spatially higher order differencing is to replace the argu¬ 
ments u,'41 and u, by u R , , and u L , , where 

«f+i = «j+i “ \ [(! - 0 ) A y+| u + (! + ^) A i+i u ] ( 3 - 3 ) 

u j+i = Uj + - |(1 - 0)Ay_iti + (1 + 0)A J + iuj . (3.4) 

Here the order of accuracy in space is determined by the value of 6. For example, if 9 = — 1, 

the scheme is fully upwind and, if 0 = 0, it is the Fromm’s scheme. When 9 = 1/2, the scheme 

is third-order and, when 9 = 1, it is the regular three-point central difference scheme. 

Various “slope” limiters are use to eliminate unwanted oscillations. A popular one is the 
“minmod” limiter; it modifies the upwind-biased interpolation (3.3)-(3.4) as follows: 
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where 


«7 + i = U 3+l 


u i+i - u ^'+ 1 + I 


(1 - 0)A J+ |iz + (1 + 0)A J+ iu 
(1 - 6 ) A i _iu+ (1 + 0)A J+ iit 


(3.5) 

(3.6) 


Ay+i« = minmod (A J+ iu,/?A 3 „iu) 
A J+ iU = minmod ^A J+ iu,/?A J+ |U^ 


with 


and 


ft 


3-0 


(3.7) 

(3.8) 


:>d(x,/?y) = sgn(x) • max jo, min [|x|,/?ysgn(x)] j (3.9) 


(3.10) 


Therefore, a higher order scheme can be obtained by simply redefining the arguments of h }±i 
in equation (3.1); i.e. 


l l +1 = u" - A 


M u f+i> u f+i) - H u f- i> u f_i) 


(3.11) 


Applied to Roe’s scheme, the second-order numerical flux in (3.11) denoted by , is 


lU _ . i,R L 

h 3 + i ~ M J+i* J+1 


/(«f + l) + /(«^l) 


/<“?+*)-/(«?+*) 


(ll* 

v J + 2 






(3.12) 


Anderson et al. [19] have demonstrated that the use of limiters is not necessary via the 
MUSCL differencing for some steady-state airfoil calculations. In reference [19], they use the 
fully one-sided unlimited second-order upwind scheme (0 = -1 in (3.4)-(3.5)) via the flux vector 
splitting of both van Leer [7] and Steger and Warming [26]. Very accurate steady-state numerical 
solutions were obtained for some typical airfoil calculations with shocks. 

For the fully upwind scheme with 0 = —1 and without limiter, uft x and u L A in (3.12) can 

Jt 2 0 • 2 

be defined as 
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or the limited version 


(3.13) 

(3.14) 


(3.15) 

(3.16) 

For the unlimited case (3.13)-(3.14), less computation is required for (3.11) than Harten’s scheme 
(2.2). However, for the limited case (3.15)-(3.16) or (3.5)-(3.6) with 9 ^ —1, an extra limiter 
is required over the Harten scheme. This is due to the asymmetric nature of the scheme. In 
general (3.12) requires one more flux evaluation than Harten scheme requires. 

IV. Osher and Chakravarthy TVD Schemes 

Instead of using the MUSCL differencing, Osher and Chakravarthy [10] started with a one 
parameter family of semi-discrete schemes with numerical flux 


ttf+i = «i +1 - rA i+f tt 


T A 

u . i = u,- 4—A, i u 
1+i 3 O J— 2 


«f + i = «y+i- §A,. + | tt 


uf, i = u 7 + - A, xu. 

3+2 3 •? 3 2 



where h J + x = h(uj, tty+i) is the numerical flux of some first-order upwind TVD scheme. For 
comparison purposes, define h J + x as in (3.2). This five-point scheme again can be second or 
third-order depending on the value of the parameter 9. It is defined the same way as before. 
The superscript + or — denotes the flux difference across the wave with positive or negative 
wave speed. To obtain a TVD scheme, Osher and Chakravarthy used the flux limiter concept 
by simply modifying the last four terms on right hand side of (4.1); i.e., their numerical flux 
denoted by h 0 f k is 

3-r 2 



The “ ~ ” shown over the flux difference values denotes the flux limited value and they are 
computed as follows 


Ay + s/ = minmod 


A y+§/ >^ A i+i/ 


(4.3a) 
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(4.3b) 


A i+if~ = minmod [ A i+i/ ,^ A J+ |/ ] 

Ay+|/ + = minmod [A y+ i/ + ,/?A y _ x/ + ] (4.3c) 

Ah/ + = minmod [A 3 _i/+,/?A y+ i/+] (4.3d) 

where the operator “minmod” is the same as (3.9) and is defined almost the same way as in 
(2.4d). The value /? is the same as in (3.10). Due to the design principle of this scheme, the 
numerical flux (4.2) requires more computation than (2.4), (2.8) or (3.12). For 0 — —1, (4.2) 
together with (4.3) requires one more limiter than Harten’s method. For 9^—1, two more 
limiters than van Leer’s method and three more limiters than Harten’s method are required for 
(4.2)-(4.3). The extra computations will become even more apparent as we extend these schemes 
to system cases. 


V. Generalization to Hyperbolic System of Conservation Laws 

The procedure considered here for generalization to system cases is to extend the scalar scheme 
to system cases so that the resulting scheme is TVD for the “locally frozen” constant coefficient 
system. To do that, we define at each point a “local” system of characteristic fields, and then 
apply the scheme to each of the m scalar characteristic equations. Here m is the dimension of the 
hyperbolic system. The extension technique is a somewhat generalized version of the procedure 
suggested by Huang [11] and Roe [12]. For simplicity of presentation, only conservation laws in 
Cartesian coordinates with equal spatial step size will be discussed. 


Consider a nonlinear hyperbolic system of conservation laws 

(5.1) 

Here U and F(U) are column vectors of m components. Let the eigenvalues of A be 
(a 1 , a 2 , ...,a m ). Denote R as the matrix whose columns are eigenvectors of A, and jf? _1 as 
the inverse of R. Define the characteristic variable of W with respect to the state U by 


d_U dF{U) 
dt dx 


A(U) 


dF(U) 
dU ' 


W = R~ 1 U. (5.2) 

In the constant coefficient case (5.1) decouples into m scalar equations for the characteristic 
variables 


dw l ,dw l 

~dT + a a^ = 0) 


(5.3) 


with ^constant and w l the elements of W. This offers a natural way of extending a scalar 
scheme to a constant coefficient system by applying it to each of the m scalar characteristic 
equations (5.3). 
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Let £/j + 1 denote some symmetric average of Uj and U J+ i such as the arithmetic mean average 
of Huang or the Roe’s average [11,12] for gas dynamics. Let a l . + i( RJ * A denote the 

quantities a 1 , R, f?~ x evaluated at £/„•, 1 . Let be the vector elements of W and let a 1 , , = 

J r 2 J+J 

~ be the component of A J+ if7 = t/y+j — Uj for the /th characteristic; i.e., 

= (5 ' 4a) 

and 

A i + h U = R i+h a i+h' ( 54b ) 

With the above notation, one can apply the scalar scheme to each of the locally defined 
characteristic variables of (5.3). Transforming back to the original variable locally, one gets 

= (5.5) 

Numerical Flux for Harten Scheme 

The corresponding numerical flux of (2.4) can be written as 

= i[f, + fV + i + tf )+4 *, +4 ]' (56a) 

where the elements of the 3> J+ i denoted by + = l,...,m are 


tin = +9 > +i ~ + K+i 

Si = minmod(u' + ia' + i,rr'_ia'_ A ) 


(5.6b) 

(5.6c) 


with xp(z) defined in (2.4b); a 1 L is (2.4f) with z = a 1 A , and 

J • 2 3 1 2 

j + 
** , 


'4 + 


(^+i - 


t+5 


1^0 

2 

0 , 


(5.6d) 


where a*. + A are the elements of (5.4a). 


Extension of the scalar TVD scheme to nonlinear system cases is not unique. Take, for 
example, the case where the numerical flux h J+ i in (2.4) is identical to h J+ i in (2.8). The 
corresponding numerical fluxes for the system case have a different form depending on which of 
the scalars hj + x one started with. If one started with (2.8), the system can be expressed as 


U? 


+1 _ 


U? 


j {> - ft”. - V )- j {‘+*U) ft" - F-i) I”*) 


where 
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and one way to define 7j is 


A 3+$- R j+t*i+l R iH 

(5.7b) 

A i+ i = diag [sg n (o‘ +i )] , 

(5.7c) 

7-j — Fj + RjGj 

(5.7d) 


with Gj = (</*, ...,g™) T and diag(z 1 ) denoting a diagonal matrix with diagonal elements z l . The 


numerical flux Hj±x can be of the form 


— ~ (l — A j+i'} (^j'+l — ?j) + 

(5.8a) 

Hj- \ ~ ?i~ 2 (* + A j- § ) ~ 0 • 

(5.8b) 


Eventhough Rj+x is needed in both (5.6) and (5.8), scheme (5.6) is preferred over (5.8) since 
Rj is not required in (5.6). In other words, for scalar cases (2.4) and (2.8) are identical, but in 
system cases (5.8a) requires more work than (5.6a). 


Numerical Flux for van Leer Differencing 
The numerical flux in this case is 

S iH = + ( 5 ' 9a > 

where Rj + x is the eigenvector of A evaluated at some symmetric average of U R + L and Uf? +L \ 
i.e., 


R,+x = R{Uf H ,Uf H ) (5.9b) 

and the elements of , x are 

^ + i = |S‘ + i |S' +i (5.9c) 

where again a 1 A and a 1 . , are evaluated at some symmetric average of U R , and U L l and 

J~v 2 3\ 2 3' 2 3 "i 2 

a; + ,=«'(c” k ,u} H ) (5.9d) 

^ .. 


s J+ 4 = «,1i (vf+i-vfm 


(5.9e) 


with 
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j = Vj+ 1 - j [(» - «)A J+ ,t/ + (1 + *)A, + Jtr] (5.9f) 

f'Z+i = C> + \ [(1 - »)A~Jt < + (1 + »)A^lf] (5.9g) 

where each element of the quantities are defined in the same way as in (3.7) and (3.8). 
Numerical Flux for Osher and Chakravarthy Schemes 
Their numerical flux for a system can be written as 


fjOC _ TT _ If A- 


C 1 + Q) R A - 

4 R i+t*j+t 


+ (L^r A+ +{LJ 1 r A + 

4 J+2 J+5 4 J 2 J-2 


Here the H„,i is the first-order numerical flux 


H i+h = 2 Fj + Fj+1 ~ R >+k 9 i+i ’ 


where the elements of + a are <f> J + i = |a*. + A |cr + i 
The elements of the A ± ’s are 

^J+ 2 ~ a j+2 a j’+2 
/.+ = a'+ o'jl 

with a'.± A = 0.5(aJ. + A ± |aJ. + A |), and a l . + A is (5.4), and 

A l - + , = minmod [a^JA 1 ^,] 

A T + i= “inmod 


A l + + , = minmod [/+ a ,/M£ a ] 
minmod [^,/Mj+j] • 


Collecting terms of f? 3+ i, (5.10a) can be rewritten as 


(5.10a) 

(5.10b) 

(5.10c) 

(5.10d) 

(5.10e) 

(5.10f) 

(5.10g) 

(5.10h) 
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(5.11) 


H°° 


H 


Fj + Fj +1 _ # 
(1-0) 




0 + 




(1 + 0 ) 


*u 


+ 4 


The terms that are common to the numerical fluxes of Harten’s H h l (5.6), van Leer’s HY*\ 

^ 3*2 3 I 2 

(5.9) and Osher and Chakravarthy’s (5.11) are Fj+i, Fj, Rj+i and a J+ x. Both H h and 

H vl require one matrix-vector multiplication + or f? J+ i$y + i. Comparing (5.6b)- 

(5.6d) and (5.9c)-(5.9g), the computational effort between H h and H VL is very competitive. 
The first three terms on the right hand side of H°9 in (5.11) require similar operation count 
as H h and H VL . The last two terms of (5.11) consist of approximately 2/5 of the computation 
of the numerical flux \ . They are the extra computation that are not present in H h and 

H vl . Therefore, for 0 7 ^ —1, H oc is approximately 40-45% more expensive to use than H h or 
H VL . For 0 = — 1, H oc is approximately 20-25% more expensive than the other two schemes. 


Concluding Remarks 

The current note addresses an aspect of implementating three existing spatially five-point 
upwind schemes for systems of hyperbolic conservation laws. The reason for choosing the above 
methods is that they represent typical design principles often used in constructing second-order 
(or higher-order) upwind TVD schemes. The current conclusion on the relative computational 
effort among the three methods mainly addresses unsteady calculations. 

It is observed that the five-point upwind TVD scheme of Osher and Chakravarthy requires 
more operations than Harten or van Leer. The reason lies mainly in the different way of arriving 
the spatially second and third-order TVD schemes. The indication is that one can achieve the 
same spatial order of accuracy as Osher and Chakravarthy by using the lower operation count 
method of van Leer or Harten. Moreover, for the forward Euler time differencing, only Harten’s 
method is second-order in time and space; the other two are first-order in time but second-order 
in space. To achieve second-order in time as well as space, additional work is required for the 
van Leer and for the Osher and Chakravarthy methods. 

For steady-state applications, aside from accuracy requirements, convergence rate of a scheme 
is equally or more important than the CPU time per time step. The performance of an individual 
scheme in terms of convergence rate is highly problem dependent. Therefore the discussion 
here requires extensive numerical testing before a clearer picture can be drawn for steady-state 
applications. 
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